%run this after full analysis
save_results = 0;
figname_out = 'other_data_summary.fig';
fname_mat_out = 'other_data_summary.mat';
fname_out = 'other_data_summary.txt';
data_folder_prefix = 'data_*';

% StartDir = '\\128.112.86.75\lithium\Publications\Charge-Diffusion\Analysis\ones_data';
% folders = dir(fullfile(StartDir, ));
% data_folder_prefix = 'data_*';
[main_dir, ~, ~] = fileparts( mfilename('fullpath') );
dir_form = fullfile(main_dir, data_folder_prefix);
folders = dir(dir_form);
folders = folders([folders.isdir]);
% folders = folders(3:end);

ds = [];

%TODO need to load all folder data to get this to work properly
for ii = 1:length(folders)
    FullPath = fullfile(folders(ii).folder, folders(ii).name);
%     if ~isempty(dir(fullfile(FullPath, '*_struct.mat')))
        dstemp = DiffuseSetSummary();
        dstemp.initialize(FullPath, [], [], 0,...
            FullPath, '', 1);
        dstemp.DiffuseSets = {}; %don't need these here...
        ds = cat(1,ds,dstemp);
%     end
    close all;
end

 %%
% first set up figure. No plotting yet. 
FigHandle_Params = figure('name', 'Amplitudes, etc');
NRows = 3;
NCols = 3;
ax_stack = cell(1, 9);

ax_stack{1} = subplot(NRows, NCols, 1);
hold on;
xlabel('T (t)');
ylabel('<n>');
grid on;
ylim([0,1]);
title('Density');

ax_stack{2} = subplot(NRows,NCols,2);
hold on;
xlabel('T (t)');
ylabel('<n>');
grid on;
ylim([0,1]);
title('All Densities');

ax_stack{3} = subplot(NRows,NCols,3);
hold on;
xlabel('T (t)');
ylabel('Grads (Dens/ 10 Latt Site)');
grid on;
title('Grads One Dir');

ax_stack{4} = subplot(NRows,NCols,4);
hold on;
xlabel('T (t)');
ylabel('Grads (Dens/ 10 Latt Site)');
grid on;
title('Grads Other Dir');

try
    ax_stack{5} = subplot(NRows, NCols, 5);
    hold on;
    xlabel('T (t)');
    ylabel('Max Amp <n^s>');
    grid on;
    ylim([0,0.2]);
    title('Max Amps, <n^s>')
catch
end

ax_stack{6} = subplot(NRows,NCols,6);
% errorbar(AllSine_Temps,OffOverAmp,OffOverAmpUnc,OffOverAmpUnc,AllSine_TempsUnc,AllSine_TempsUnc,'.');
hold on;
% plot([0,max(Ts)],[0.8,0.8],'k--');
xlabel('T (t)');
ylabel('Offset/Amp');
grid on;
ylim([-1,1]);
title('Offsets/Amp ~= Non-linearity, <n_1>')

ax_stack{7} = subplot(NRows,NCols,7);
% errorbar(AllSine_Temps,OffOverAmp,OffOverAmpUnc,OffOverAmpUnc,AllSine_TempsUnc,AllSine_TempsUnc,'.');
hold on;
% plot([0,max(Ts)],[0.8,0.8],'k--');
xlabel('T (t)');
ylabel('Offset');
grid on;
ylim([-0.1,0.1]);
title('Offsets, <n_1>')

%TODO: add these to analysis
ax_stack{8} = subplot(NRows, NCols, 8);
hold on;
xlabel('T (t)');
ylabel('Gamma (KHz)');
grid on;
% ylim([-0.1,0.1]);
title('Gamma')

ax_stack{9} = subplot(NRows, NCols, 9);
hold on;
xlabel('T (t)');
ylabel('D (latt^2/ms)');
grid on;
% ylim([-0.1,0.1]);
title('D')


suptitle('Parameters Summary');

cm = get_color_map(length(ds), [], [0, 1, 1], []);

% some numbers for all sets
fidelity = 0.97;
fidelity_unc = 0.02;

% sort sets by temperature
Ts = [ds.Temp];
[Ts,I] = sort(Ts);
for ii = 1:length(ds)
    
    % extract data from Diffuse sets
    ds_temp = ds(I(ii));
    Ts = ds_temp.Temp;
    TsUnc = ds_temp.TempSD;
    Dens = ds_temp.Dens;
    DensUnc = ds_temp.DensSD;
    
    AllDens_Temps = ds_temp.AllTemps;
    AllDens_TempsUnc = ds_temp.AllTempsUnc;
    
    AllDens = ds_temp.AllDens;
    AllDensUnc = ds_temp.AllDensUnc;
    
    dens_corr = mean(AllDens) / fidelity;
    dens_corr_unc = dens_corr * sqrt((fidelity_unc / fidelity)^2 + (std(AllDens)/ mean(AllDens) )^2);
    AllGrads = ds_temp.AllDensGrads;
    AllGradsUnc = ds_temp.AllDensGradsUnc;
    AllAmps = ds_temp.AmpsHydro;
    AllAmpsUnc = ds_temp.AmpsStdErrsHydro;

    % Data from individual DataFolders
    % All sine fits, for each individual DataFolder
    AllSineFpCell = ds_temp.AllSineFitParams;
    AllSineErrCell = ds_temp.AllSineFitErrs;

    AllSineFp = [];
    AllSineErr = [];
    AllSine_Temps = [];
    AllSine_TempsUnc = [];
    for jj = 1:length(AllSineFpCell)
        AllSine_Temps = cat( 2,AllSine_Temps,...
            AllDens_Temps(jj) * ones(1, length( AllSineFpCell{jj} )) );
        AllSine_TempsUnc = cat(2, AllSine_TempsUnc,...
            AllDens_TempsUnc(jj) * ones(1, length(AllSineFpCell{jj})));
        AllSineFp = cat(2, AllSineFp, AllSineFpCell{jj});
        AllSineErr = cat(2, AllSineErr, AllSineErrCell{jj});
    end

    SineFp_FirstTime = zeros(4, length(AllSineFpCell));
    SineErr_FirstTime = zeros(4, length(AllSineFpCell));
    for jj = 1:length(AllSineFpCell)
        SineFp_FirstTime(:,jj) = AllSineFpCell{jj}(:,1);
        SineErr_FirstTime(:,jj) = AllSineErrCell{jj}(:,1);
    end

    OffOverAmp = SineFp_FirstTime(4,:) ./ abs(SineFp_FirstTime(2,:));
    OffOverAmpUnc = abs(OffOverAmp) .* sqrt((SineErr_FirstTime(2,:) ./ SineFp_FirstTime(2,:)).^2 + ...
        (SineErr_FirstTime(4,:) ./ SineFp_FirstTime(4,:)).^2);

    Offset = SineFp_FirstTime(4,:);
    OffsetUnc = SineErr_FirstTime(4,:);
    
    
    
    % plot results
    errorbar(ax_stack{1}, Ts, Dens, DensUnc, DensUnc, TsUnc, TsUnc,...
        '.', 'Color', cm(ii, :));
    errorbar(ax_stack{2}, AllDens_Temps, AllDens, AllDensUnc, AllDensUnc,...
        AllDens_TempsUnc, AllDens_TempsUnc,...
        '.', 'color', cm(ii, :));
    errorbar(ax_stack{3}, AllDens_Temps, AllGrads(1,:)*10, AllGradsUnc(1,:)*10,...
        AllGradsUnc(1,:)*10, AllDens_TempsUnc,AllDens_TempsUnc,...
        '.', 'color', cm(ii, :));
    errorbar(ax_stack{4}, AllDens_Temps, AllGrads(2,:) * 10,AllGradsUnc(2,:) * 10,...
        AllGradsUnc(2,:) * 10, AllDens_TempsUnc, AllDens_TempsUnc,...
        '.', 'color', cm(ii, :));
    errorbar(ax_stack{5}, AllDens_Temps, AllAmps, AllAmpsUnc, AllAmpsUnc,...
            AllDens_TempsUnc, AllDens_TempsUnc, '.', 'color', cm(ii, :));
    errorbar(ax_stack{6}, AllDens_Temps,OffOverAmp,OffOverAmpUnc,...
        '.', 'color', cm(ii, :));
    errorbar(ax_stack{7}, AllDens_Temps,Offset,OffsetUnc,...
        '.', 'color', cm(ii, :));
end

% plot desired density lines
plot(ax_stack{1}, [0, max(Ts)], [0.8,0.8],'k--');
plot(ax_stack{2}, [0,max(Ts)],[0.8,0.8],'k--');

% get parameters for everything
%plot results 
%overall densities and temperatures for each collection of DiffuseSets
Ts = [ds.Temp];
[Ts,I] = sort(Ts);
TsUnc = [ds(I).TempSD];
Dens = [ds(I).Dens];
DensUnc = [ds(I).DensSD];

%densities, temperatures for each DiffuseSet, etc.
AllDens_Temps = [ds(I).AllTemps];
AllDens_TempsUnc = [ds(I).AllTempsUnc];

AllDens = [ds(I).AllDens];
AllDensUnc = [ds(I).AllDensUnc];
fprintf('Density = %0.2f +/- %0.2f\n', mean(AllDens), std(AllDens));

dens_corr = mean(AllDens) / fidelity;
dens_corr_unc = dens_corr * sqrt((fidelity_unc / fidelity)^2 + (std(AllDens)/ mean(AllDens) )^2);
fprintf('Fidelity Estimate = %0.2f +/- %0.2f\n', fidelity, fidelity_unc);
fprintf('Corrected Density = %0.2f +/- %0.2f\n', dens_corr, dens_corr_unc);
AllGrads = [ds(I).AllDensGrads];
fprintf('Grad 1 = %0.2f +/- %0.2f\n', mean(AllGrads(1,:)), std(AllGrads(1,:)));
fprintf('Grad 2 = %0.2f +/- %0.2f\n', mean(AllGrads(2,:)), std(AllGrads(2,:)));
AllGradsUnc = [ds(I).AllDensGradsUnc];
AllAmps = [ds(I).AmpsHydro];
AllAmpsUnc = [ds(I).AmpsStdErrsHydro];
fprintf('Amp = %0.2f +/- %0.2f\n', mean(AllAmps), std(AllAmps));

%All sine fits, for each individual DataFolder
AllSineFpCell = [ds(I).AllSineFitParams];
AllSineErrCell = [ds(I).AllSineFitErrs];

AllSineFp = [];
AllSineErr = [];
AllSine_Temps = [];
AllSine_TempsUnc = [];
for ii = 1:length(AllSineFpCell)
    AllSine_Temps = cat(2,AllSine_Temps,AllDens_Temps(ii)*ones(1,length(AllSineFpCell{ii})));
    AllSine_TempsUnc = cat(2,AllSine_TempsUnc,AllDens_TempsUnc(ii)*ones(1,length(AllSineFpCell{ii})));
    AllSineFp = cat(2,AllSineFp,AllSineFpCell{ii});
    AllSineErr = cat(2,AllSineErr,AllSineErrCell{ii});
end

SineFp_FirstTime = zeros(4,length(AllSineFpCell));
SineErr_FirstTime = zeros(4,length(AllSineFpCell));
for ii = 1:length(AllSineFpCell)
    SineFp_FirstTime(:,ii) = AllSineFpCell{ii}(:,1);
    SineErr_FirstTime(:,ii) = AllSineErrCell{ii}(:,1);
end

OffOverAmp = SineFp_FirstTime(4,:) ./ abs(SineFp_FirstTime(2,:));
OffOverAmpUnc = abs(OffOverAmp) .* sqrt((SineErr_FirstTime(2,:)./SineFp_FirstTime(2,:)).^2 + ...
    (SineErr_FirstTime(4,:) ./ SineFp_FirstTime(4,:)).^2);

Offset = SineFp_FirstTime(4,:);
OffsetUnc = SineErr_FirstTime(4,:);



if save_results
    savefig(FigHandle_Params,figname_out);
    save(fname_mat_out, 'Ts', 'TsUnc', 'Dens', 'DensUnc', 'AllDens', 'AllDensUnc',...
        'fidelity', 'fidelity_unc', 'dens_corr', 'dens_corr_unc',...
        'AllGrads', 'AllGradsUnc', 'AllAmps', 'AllAmpsUnc',...
        'AllSineFp', 'AllSineErr', 'AllSine_Temps', 'AllSine_TempsUnc',...
        'SineFp_FirstTime', 'SineErr_FirstTime');
end